Scaling behavior of the disordered contact process 
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The one-dimensional contact process with weak to intermediate quenched disorder in its trans- 
mission rates is investigated via quasi-stationary Monte Carlo simulation. We address the contested 
questions of both the nature of dynamical scaling, conventional or activated, as well as of univer- 
sality of critical exponents by employing a scaling analysis of the distribution of lifetimes and the 
quasi-stationary density of infection. We find activated scaling to be the appropriate description for 
all disorder strengths considered. Critical exponents are disorder dependent and approach the values 
expected for the limit of strong disorder as predicted by strong-disorder renormalization group anal- 
ysis of the process. However, even for the strongest disorder under consideration no strong-disorder 
exponents are found. 



43 

o 

B 

i 



o 
o 



> 
q 

On 
O 
00 

o 



X 



The critical behavior of systems with quenched ran- 
domness has been the subject of interest for a long time. 
Over the past decades, their investigation has revealed 
rich behavior including the existence of new phases [l|, 
novel fixed points 0], and unconventional scaling prop- 
erties 0. 

Recently, attention has turned towards the influence 
of disorder on stochastic many-particle systems with a 
phase transition into an absorbing state owing to their 
relevance in physics, chemistry and biology [4]. In partic- 
ular, the contact process (CP) p|, a paradigmatic model 
for the stochastic spreading of an infectious disease, has 
been investigated as a representative for the prominent 
universality class of directed percolation (DP). Interest 
in the influence of disorder on this process was sparked 
by the surprising lack of experimental observation of DP 
behavior in real systems, for which disorder may be re- 
sponsible Q . With a recent study presenting convincing 
evidence of DP critical behavior in turbulent liquid crys- 
tals 0] , an understanding of the effects of disorder is more 
relevant than ever. 

Initial Monte Carlo (MC) studies of the disordered 
CP (DCP) found continuously varying dynamical criti- 
cal exponents assuming conventional scaling 0, [{J [To| . 
Recently, deep insight was gained through a strong- 
disorder renormalization group study of the DCP which 
revealed an infinite-randomness fixed point (IRFP) for 
sufficiently strong disorder in close analogy to the ran- 
dom transverse-field Ising model [Til ] . While this makes 
the strong-disorder limit of the process well understood 
and predicts new strong-disorder exponents as well as 
an unconventional "activated" form of dynamical scal- 
ing, the behavior in the weak- and intermediate-disorder 
regime remains a subject of debate fl2l [l3j . Initial 
density-matrix renormalization group (DMRG) studies 
were not able to convincingly distinguish the two al- 
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tentative dynamical scaling scenarios, conventional or 
activated, and reported critical exponents continuously 
changing with disorder strength [12j. In contrast, a re- 
cent MC study reported evidence for activated scaling 
with strong-disorder exponents for all disorder strengths 

Static scaling in the DCP was found to be of conven- 
tional form [3] as predicted by the strong-disorder renor- 
malization study [ll|. However, there exists conflicting 
evidence as to the universality of the exponents for weak 
and intermediate disorder with the literature repor ting 
both disorder-dependent [TJ, [H| and strong-disorder |l3j ] 
exponents. 

In this paper, we aim to address both the question of 
the type of dynamical scaling as well as of universality of 
exponents in the weak- and intermediate-disorder regime 
by considering the scaling of the distribution of lifetimes, 
P(t), of the Id DCP obtained from quasi-stationary MC 
simulation. This is motivated by the fact that an analysis 
of the scaling behavior of entire distributions promises to 
yiel d clearer results as compared to the scaling of means 
[I3.ll6l|. Further, in dynamic single-seed MC simulations 
employed for the DCP in the past [TH , the question 
of whether the long-time limit of the process had been 
reached was frequently contested. In contrast, quasi- 
stationary simulations offer a clear means of ensuring 
this: a true stationary average whose convergence can 
be monitored. 

In the clean CP (without disorder) defined on a lattice, 
sites represent the individuals of a population which can 
be in two possible states, susceptible or infected. An 
infected site attempts to spread its infection to nearest 
neighbors at rate A, while recovery is spontaneous at rate 
e = 1. In the limit t — > oo for an infinite system, there 
exist two distinct states: an active one where a finite den- 
sity p of infected sites remains and a non-active regime in 
which the system ultimately gets trapped in an absorb- 
ing state with no infected sites remaining. The system 
undergoes a continuous phase transition between these 
two phases at a critical rate A — A° = A = with order 
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parameter p 0,[T3|- 

Starting from a fully-infected system, the density of 
infected sites initially relaxes while spatial correlations 
grow towards the size of the system and temporal corre- 
lations decay. Once the correlation length becomes com- 
parable to this size, the process enters a quasi-stationary 
(QS) regime. This metastable state is characterized by 
a non-zero time-independent transition rate to the ab- 
sorbing state. Given that no true non-trivial stationary 
state can exist in a finite system, in these cases QS av- 
erages are commonly used as a proxy in the CP and al- 
lied models. Ultimately, the clean CP is bound to enter 
the absorbing state, the approach of which is character- 
ized by an exponentially decaying probability of survival, 
P s {t) ~ exp(— t/r), with a characteristic lifetime r. In 
the active state, A > 0, this lifetime is known to obey 
finite-size scaling, 

r - N z exp (N d A"- Ld ) , (1) 

where z and v± are critical exponents while N denotes 
the size of the system. Further, the QS density of infected 
sites, p, obeys a similar scaling form characterized by the 
universal scaling function F, i.e. 

p~N- x F(AN- 1/v± } , (2) 

with x = f3/v±, where f3 is the order parameter expo- 
nent in the infinite system defined via the behavior of 
the control parameter in the vicinity of the critical point, 

p ~ a 13 0, 03 . 

Turning to the disordered process, we follow Ref. [l3[ 
and incorporate quenched randomness into the transmis- 
sion rates of individual sites i, Aj, which are drawn from 
the bimodal distribution 

P(Ai) = (1 - P )S(\ l - A) + pS(Xi - cA) , (3) 

where p controls the concentration of impurities while 
< c < 1 characterizes the strength of the disorder. As 
c < 1, a particular realization of the disorder contains 
a concentration p of randomly arranged impurity sites 
which are less active than the surrounding sea of host 
sites. As a consequence, there now exists a new dirty 
critical point at a rate A c > A". Also, observables may 
in general take different values between different realiza- 
tions leading to disorder-induced distributions such as 
P(t) for the lifetime. 

Scaling predictions for P(t) exist and differ for the two 
alternative scaling scenarios, i.e. conventional and acti- 
vated scaling. In the former case, the relevant variable is 
r and its average over disorder, (r), is expected to obey 
a scaling form analogous to Eq. {]} albeit with a possibly 
disorder-dependent dynamical exponent z. Accordingly, 
the appropriate scale-invariant combination of variables 
is tN~ z and the lifetime distribution at criticality is ex- 
pected to scale as 

P{t) = N- z P(tN~ z ) , (4) 



where P is a universal scaling function. 

Systems that exhibit activated scaling are charac- 
terized by a strong dynamical anisotropy: the typical 
length-scale is related to the logarithm of the typical 
timescale, thus rendering the dynamical exponent for- 
mally infinity. This reflects the notion that the very 
broad distributions for observables are better described 
by their geometric rather than arithmetic means (2j. 
For the lifetime, this leads to a scaling combination 
A~* ln(r) and a corresponding scaling relation [Til ]. 

P(ln(r)) ~ AT"*P(iV-* ln(r)) , (5) 

where "J is an activated scaling exponent [cf. Eq. Q]. 
For an IRFP, which is known to control the critical be- 
havior of the DCP for sufficiently strong disorder, the ex- 
ponent takes the value \P = i . This type of fixed point is 
characterized by an extreme dynamical anisotropy, ultra- 
slow dynamics and distributions whose width diverges 
with system size 0, [HI . 

The QS density p is expected to follow a conventional 
scaling form analogous to Eq. ([2]) with an exponent x 
which differs from the clean DP value and, for sufficiently 
strong disorder, is predicted to be x = 3 ~/^ « 0.19 [111 ]. 

In order to check the above scaling relations for the 
DCP, the QS state can be investigated numerically. Anal- 
ysis of this metastable state in computer simulations has 
proved to be notoriously difficult in the past. Commonly, 
the time-dependent density of infected sites conditioned 
on survival, which becomes stationary in the QS regime, 
is investigated [l7| . Problematically though, it is neither 
clear a priori at what time an observable like the average 
density of infected sites p has converged to its QS value 
nor when the QS state starts to decay due to finite size 
effects [lj| . Therefore, a range of alternative approaches 
have been proposed which enable an observation of this 
metastable regime (see Ref. [2(| and references therein). 
Here, we employ the QS simulation method [2(| which al- 
lows a direct sampling of the QS state by eliminating the 
absorbing state and redistributing its probability mass 
over the active states according to the history of the pro- 
cess. The modified process possesses a true stationary 
state which corresponds to the original QS state and al- 
lows a precise measurement of QS observables. Generally, 
the method has proved to be efficient with fast and reli- 
able convergence after optimization of history sampling 
parameters. As demonstrated in Ref. [2l[, the lifetime of 
the QS state can be determined as the inverse of the rate 
of attempts by the system to enter the absorbing state. 

In order to investigate the validity of the two differ- 
ent scaling scenarios for the lifetime distribution at the 
dirty critical point, simulation data for disorder strengths 
c = 0.2,0.4,0.6 and 0.8 at a concentration p = 0.3 were 
considered at the critical rates reported in Ref. [13j . Sys- 
tem sizes of N — 16, 32, 64 and 128 sites were investigated 
with data sampled from no less than 10 4 disorder real- 
izations per system size and QS simulation times of up 
to 10 9 time steps. Given that probability density func- 
tions (PDFs) are difficult to obtain from simulation data 
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(a) c = 0.2, * = 0.29 and z = 3.44 




0.0 0.5 1.0 1.5 2.0 2.5 3.0 

\n(tL~ z ) 

(b) c = 0.8, * = 0.22 and z = 1.65 

FIG. 1: (Color online) Scaling collapse for the distributions of 
lifetimes r for system sizes N = 16 (o), 32 (□), 64 (O), 128 (V) 
for disorder strengths c = 0.2 (a) and c = 0.8 (b) according 
to both activated (main panel) and conventional (inset) scal- 
ing predictions. Scaling exponents were determined from the 
finite-size scaling of the means (ln(r)) or (r), respectively, as 
described in the text. 



(as binning procedures have to be used which may intro- 
duce artifacts), we perform scaling on cumulative den- 
sity functions (CDFs), F T (t) = f*P(t')dt'. The scaling 
properties of these can be derived by starting from the 



two scaling forms for conventional and activated scaling, 
given by Eqs. (|3]) and ([5]), respectively, where for the for- 
mer case 



F T (t) 



N~ z P{t'N- z ) dt' = Q{tN~ z ) 



(6) 



with Q(x) a new scaling function. An analogous expres- 
sion follows for the case of activated scaling with r re- 
placed by ln(r). 

As displayed in Fig. [1] the resulting CDFs were col- 
lapsed onto each other according to the two possible scal- 
ing scenarios (main panels). In order to achieve a fair 
comparison, logarithmic variables were used for the con- 
ventional scaling case [l6[ . As illustrated in Fig. [21 the 
exponents '3/ and z were determined from a power-law 
fit to the appropriate scaling forms for the mean in the 
two scenarios, i.e. (ln(r)) ~ TV* and (r) ~ N z . Insets 
show the alternative collapse using PDFs as discussed 
above which requires the use of histograms. There, the 
size of bins was chosen in order to minimize noise and 
the smooth curve was obtained by Gaussian broadening 
of individual data points. 

For the case of strongest disorder (c = 0.2, A c = 5.24), 
least-squares fitting gave exponents ^ c =o.2 = 0.29(2) and 
z c =o.2 = 3.44(3) for the two scaling scenarios, respec- 
tively. Data collapses for the distributions are shown 
in Fig. 1(a) for both activated (top panel) and conven- 
tional scaling (bottom panel). From these results, we 
judge that the activated scaling scenario provides a bet- 
ter fit to the data. In particular, while the collapse is not 
perfect, it is not found to exhibit any systematic trends 
which would hint at a fundamental inconsistency with the 
scaling form. This is also confirmed by the inset which 
shows the corresponding collapse of the PDF. Generally, 
while the fit is excellent for small to medium values of r, 
it gets worse with increasing r. We attribute this to the 
fact that large values of r correspond to rare events caus- 
ing the tail of the distribution to have been sampled at a 
comparatively poorer density than the bulk which results 
in stronger fluctuations. Considering the collapse using 



the conventional scaling form [bottom panel of Fig. 1(a) 
on the other hand produces a clear trend of shifting of 
distributions between different system sizes indicating a 
worse collapse as compared to the previous case for both 
CDF and PDF. 

Looking at an analogous analysis for the case of weak- 
est disorder (c = 0.8, A c = 3.525) in Fig. 1(b) the two 
scaling scenarios become harder to differentiate. A col- 
lapse using the measured exponents ^ c =o.s = 0.22(2) and 
z c= o.s = 1.65(3) appears to work similarly well in both 
cases but the quality of collapse is too poor to allow any 
definitive judgment. A close look reveals a systematic 
trend of shifting curves for the case of conventional scal- 
ing [bottom panel in Fig. l(b)| while crossings appear to 
be less syste matic in the activated scaling picture [top 



panel in Fig. 



Therefore, slight preference may be 



given to the activated scaling scenario. 

In addition, the intermediate cases of both c 
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FIG. 2: (Color online) Example of the scaling of the average 
(ln(Y)) (left) and (r) (right) with system size for the case of 
strongest disorder (c = 0.2) used to extract the exponents 
and z. 
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FIG. 3: (Color online) The quasi-stationary density p as a 
function of system size N at criticality for both the strongest 
(top, c = 0.2) and the weakest (bottom: c = 0.8) disorder 
strength. Dashed lines are linear least-squares fits which gave 
slopes of x c= o.2 = 0.22(1) and x c =o.8 = 0.25(1) while the 
solid line is a guide to the eye with slope ^strong = 0.19, the 
expected exponent in the strong-disorder limit. 



and c = 0.6 (not shown) were considered in an analogous 
fashion and showed an excellent activated scaling collapse 
of similar quality as for the case of c = 0.2. Collapsing 
data for all disorder strengths considered on the same 
plot shows no universality of the scaling function between 
different disorder strengths, i.e. it appears to be disorder 
dependent. 

Finite-size scaling of the QS density p(N), as shown 
in Fig. 02 is found to be conventional with a disorder- 
dependent exponent x = x(c). The conventional scal- 
ing form is in line with previous investigations [14| and 
theoretical predictions [ll|. Moreover, the fact that we 
again find continuously varying exponents gives addi- 
tional credibility to the scaling picture presented above. 
Generally, both the exponent 'J' and the exponent x are 
found to approach their values predicted from strong- 
disorder renormalization [ill ] with increasing strength of 
disorder starting from their homogeneous values as shown 
in Fig. |U For the strongest disorder under consideration 



FIG. 4: (Color online) Critical exponents (top panel) and 
x (bottom panel) as a function of disorder strength c where 
dashed lines show the values in the limit of strong disorder. 



(c = 0.2), both exponents are found to still be well away 
from their predicted strong-disorder values. 

The above findings suggest that for strong enough dis- 
order, activated scaling captures the behavior of the dis- 
ordered CP well compared to a conventional scaling pic- 
ture. For weak disorder however, no such clear conclusion 
can be made. Further, associated critical exponents ap- 
pear to change smoothly from their clean DP values ap- 
proaching their values characteristic of an IRFP asymp- 
totically with increasing disorder strength. While this 
conclusion appears to be in conflict to that presented in 
Ref. [HI , the authors of that reference do discuss doubts 
about the universality of exponents and cannot exclude 
the possibility of a change with disorder strength. Wc 
have re-analyzed some of the data presented there and 
found it to be compatible with our exponent values. 

There exist three possible explanations compatible 
with these findings. First, a continuous line of fixed 
points, one for each strength of disorder, could be present 
which for sufficiently strong disorder turns to an attrac- 
tive flow into the IRFP as suggested in Ref. [ll|. Sec- 
ond, identical and numerically indistinguishable behavior 
could be explained by a crossover between the clean DP 
fixed point and the IRFP where effective exponents are 
observed at intermediate disorder strengths due to the 
influence of both fixed points. This has been observed 
in several disordered equilibrium systems as discussed in 
e.g. Refs. [13, HU . Lastly, in principle, the observed be- 
havior could also be explained by an abrupt jump from 
the clean DP exponents to those of the IRFP obscured 
by finite-size corrections. 

The last option we feel can be excluded in light of the 
facts that perturbative series expansions (cf. Ref. [l5[) do 
not show a jump in exponents and that no evidence for 
strong corrections to finite-size scaling were observed by 
us. At the same time, the other two scenarios are compat- 
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ible with our and most other results but cannot be safely 
distinguished by numerical investigation alone without 
an established theoretical framework for the crossover in 
the DCP. 
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